* Webb, Linn and Lebo Procedure for independently-generated unit roots DGP
* Coverage

// Basics
clear all
set matsize 11000

// Set working directory

parallel setclusters 8

capture program drop NG
program define NG, rclass

// Monte Carlo Basics
version 12.1
syntax [, obs(integer 1) phi_y(real 0.5) phi_x(real 0.5)]
drop _all
set obs `obs'

// Time Series Basics
egen time = fill(0 1 2)
tsset time

* Assign bounds from WLL AJPS 2020 Table 2, alpha = 0.05, k = 1
if `obs' == 25 {
	local lb = 1.25
	local ub = 3.79
} 
else if `obs' == 50 {
	local lb = 1.09
	local ub = 3.72
}
else if `obs' == 75 { // 75 observations
	local lb = 1.06
	local ub = 3.73
}
else if `obs' == 150 { // 150 observations
	local lb = 1.01
	local ub = 3.68
}
else if `obs' == 500 { // 500 observations
	local lb = 0.99
	local ub = 3.66
}

else { // 1000 observations
	local lb = 0.99
	local ub = 3.62 
}
// Data Generating Process

* Unit root predictor
gen x = rnormal()
replace x = `phi_x' * L.x + rnormal() if time > 0

* Generate outcome
gen y = rnormal()
replace y = `phi_y' * L.y + rnormal() if time > 0

// Model Estimation


* Estimate partial adjustment 
reg y L.y x
return scalar lrm_pa = _b[x]/(1-_b[L1.y])
nlcom (1e+3)* (_b[x]/(1-_b[L1.y]))
matrix define var_lrm_pa = r(V)
scalar var_lrm_pa_se  = sqrt(var_lrm_pa[1,1])/(1e+3)
scalar lrm_pa_test = ((_b[x]/(1-_b[L1.y])) - 0) /var_lrm_pa_se
return scalar ftr_pa = abs(lrm_pa_test) <= `lb'
return scalar inc_pa = abs(lrm_pa_test) > `lb' & abs(lrm_pa_test) < `ub'
return scalar rej_pa = abs(lrm_pa_test) >= `ub'

* Estimate GECM 
reg D.y L.y D.x L.x
return scalar lrm_gecm = -_b[L1.x] / _b[L1.y]
nlcom (-_b[L1.x] / _b[L1.y])* (1e+3)
matrix define var_lrm_gecm = r(V)
scalar var_lrm_gecm_se  = sqrt(var_lrm_gecm[1,1]) / (1e+3)
scalar lrm_gecm_test = abs(((-_b[L1.x] / _b[L1.y]) - 0) /var_lrm_gecm_se)
return scalar ftr_gecm = lrm_gecm_test <= `lb'
return scalar inc_gecm = lrm_gecm_test > `lb' & lrm_gecm_test < `ub'
return scalar rej_gecm = lrm_gecm_test >= `ub'

* Estimate ADL 
reg y L.y x L.x
return scalar lrm_adl = (_b[x] + _b[L1.x]) /(1-_b[L1.y])
nlcom ((_b[x] + _b[L1.x]) /(1-_b[L1.y])) * (1e+3)
matrix define var_lrm_adl = r(V)
scalar var_lrm_adl_se  = sqrt(var_lrm_adl[1,1]) / (1e+3)
scalar lrm_adl_test = abs((((_b[x] + _b[L1.x]) /(1-_b[L1.y])) - 0) /var_lrm_adl_se)
return scalar ftr_adl = lrm_adl_test <= `lb'
return scalar inc_adl = lrm_adl_test > `lb' & lrm_adl_test < `ub'
return scalar rej_adl = lrm_adl_test >= `ub'

* Estimate ADL(2,2)
reg y L.y L2.y x L.x L2.x
return scalar lrm_adl22 = (_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y]))*(1e+3)
matrix define var_lrm_adl22 = r(V)
scalar var_lrm_adl22_se  = sqrt(var_lrm_adl22[1,1]) / (1e+3)
scalar lrm_adl22_test = abs((((_b[x] + _b[L1.x] + _b[L2.x]) /(1-_b[L1.y]-_b[L2.y])) - 0) /var_lrm_adl22_se)
return scalar ftr_adl22 = lrm_adl22_test <= `lb'
return scalar inc_adl22 = lrm_adl22_test > `lb' & lrm_adl22_test < `ub'
return scalar rej_adl22= lrm_adl22_test >= `ub'

* Estiamte ADL(3,3)
reg y L.y L2.y L3.y x L.x L2.x L3.x
return scalar lrm_adl33 = (_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])
nlcom ((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) * (1e+3)
matrix define var_lrm_adl33 = r(V)
scalar var_lrm_adl33_se  = sqrt(var_lrm_adl33[1,1]) / (1e+3)
scalar lrm_adl33_test = abs((((_b[x] + _b[L1.x] + _b[L2.x] + _b[L3.x]) /(1-_b[L1.y]-_b[L2.y]-_b[L3.y])) - 0) /var_lrm_adl33_se)
return scalar ftr_adl33 = lrm_adl33_test <= `lb'
return scalar inc_adl33 = lrm_adl33_test > `lb' & lrm_adl33_test < `ub'
return scalar rej_adl33= lrm_adl33_test >= `ub'
end

